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Abstract 

In  this  article  we  introduce  the  multi-domain  hybrid  Spectral- WENO  method  aimed 
at  the  discontinuous  solutions  of  hyperbolic  conservation  laws.  The  main  idea  is  to 
conjugate  the  non-oscillatory  properties  of  the  high  order  weighted  essentially  non- 
oscillatory  (WENO)  finite  difference  schemes  with  the  high  computational  efficiency 
and  accuracy  of  spectral  methods.  Built  in  a  multi-domain  framework,  subdomain 
adaptivity  in  space  and  time  is  used  in  order  to  maintain  the  solutions  parts  exhibit¬ 
ing  high  gradients  and  discontinuities  always  inside  WENO  subdomains  while  the 
smooth  parts  of  the  solution  are  kept  in  spectral  ones.  A  high  order  multi-resolution 
algorithm  by  Ami  Harten  is  used  to  determine  the  smoothness  of  the  solution  in 
each  subdomain.  Numerical  experiments  with  the  simulation  of  compressible  flow 
in  the  presence  of  shock  waves  are  performed. 

Key  words:  Spectral,  WENO,  Multi-resolution,  Multi-Domain,  Hybrid, 
Conservation  Laws 


1  Introduction 


The  fine  scale  and  delicate  structures  of  physical  phenomena  related  to  tur¬ 
bulence  demand  the  utilization  of  high  order  methods  when  performing  nu¬ 
merical  simulations.  Spectral  methods  are  nondispersive  and  nondissipative 
and,  therefore,  well  htted  to  this  task  when  the  solutions  involved  are  smooth. 
However,  in  the  modeling  of  compressible  turbulent  flows  by  means  of  the 
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inviscid  Euler  Equations,  the  development  of  finite  time  discontinuities  gen¬ 
erates  global  0(1)  oscillations,  known  as  the  Gibbs  phenomenon,  causes  loss 
of  accuracy  and  numerical  instability.  Filtering  of  small  scales  has  been  used 
to  stabilize  the  spectral  scheme  in  shock  calculations  [14],  however  the  Gibbs 
oscillations  still  remain  in  the  solution. 

In  hgure  1,  the  density  of  the  Lax  shock  tube  problem  with  Riemann  initial 
conditions  was  obtained  with  a  Ghebyshev  collocation  method  stabilized  with 
a  16_th  order  Exponential  hlter.  Note  that  the  oscillations  also  occur  at  a 
smaller  scale  on  the  edges  of  the  rarefaction  wave,  due  to  the  discontinuities 
at  the  derivative  of  the  density.  The  use  of  heavier  hltering  is  not  recommended 
as  it  would  remove  the  hne  scale  physical  structures  which  are  important  for 
the  overall  hdelity  of  the  simulation. 


Fig.  1.  The  density  of  the  Lax  problem  as  computed  by  (Left)  the  Ghebyshev  col¬ 
location  solution  with  the  Gibbs  oscillations  and  (Right)  the  fifth-order  character¬ 
istic-wise  WENO  finite  difference  scheme.  The  exact  solution  is  represented  by  the 
solid  line. 

Reconstruction  techniques  such  as  the  direct  and  inverse  Gegenbauer  expan¬ 
sions  (see  [4,29]  and  references  contained  therein)  have  achieved  some  success 
as  a  postprocessing  treatment  to  remove  the  Gibbs  oscillations.  These  tech¬ 
niques  followed  the  achievement  of  relevant  theoretical  results  demonstrating 
convergence  properties  of  spectral  approximations  of  discontinuous  solutions. 
For  instance,  Gottlieb  and  Tadmor  [5]  proved  that,  for  linear  problems,  the  mo¬ 
ments  of  the  numerical  solution  computed  by  spectral  methods  are  spectrally 
accurate.  Lax  [9]  had  argued  that  high  order  information  about  the  solution 
is  contained  in  high  resolution  schemes,  even  for  nonlinear  problems.  Hence 
highly  accurate  essentially  non-oscillatory  solution  can  be  extracted  from  the 
seemingly  oscillatory  noisy  data.  Tadmor  [11]  showed  the  convergence  of  spec¬ 
tral  methods  for  nonlinear  scalar  equations.  Nevertheless,  the  rapid  growth 
rate  of  the  Gegenbauer  polynomials  cause  several  numerical  problems  which 
are  difficult  to  overcome.  Very  often,  the  domain  must  be  subdivided  in  or¬ 
der  to  apply  the  Gegenbauer  reconstruction  when  dealing  with  complex  and 
localized  flow  structures. 
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These  numerical  difficulties  associated  with  the  global  spectral  methods  are 
some  of  the  main  reasons  why  local  nonlinear  adaptive  shock  capturing  hnite 
differences  have  been  the  method  of  choice  when  dealing  with  nonlinear  conser¬ 
vation  laws.  Among  these,  the  most  commonly  used  are  the  Essentially  Non- 
Oscillatory  schemes  (ENO)  [22],  In  order  to  avoid  numerical  oscillations,  ENO 
schemes  bias  the  local  stencil,  discarding  interpolations  across  discontinuities 
when  computing  the  tendencies  of  the  numerical  solution.  The  Weighted  Es¬ 
sentially  non-Oscillatory  (WENO)  method  ([24])  is  an  improvement  over  the 
ENO  method  with  the  introduction  of  a  convex  combination  of  all  the  avail¬ 
able  stencils.  WENO  achieves  optimal  order  of  accuracy  at  smooth  parts  of 
the  solution  with  the  same  stencil  size  of  ENO.  Nevertheless,  the  intrinsic 
numerical  dissipation  of  WENO  schemes,  although  necessary  to  properly  cap¬ 
ture  shock  waves,  eventually  damp  relevant  small  scales,  even  when  these  are 
smooth  components  of  the  solution.  Figure  1  also  shows  a  solution  of  the  same 
Lax  problem  obtained  with  a  hfth  order  characteristic-wise  WENO  hnite  dif¬ 
ference  scheme.  Note  the  smearing  of  the  shock,  the  contact  interface  and  of 
the  edges  of  the  rarefaction  wave.  The  numerical  dissipation  can  be  reduced 
by  increasing  the  number  of  points,  as  well  as  the  order  of  the  WENO  scheme 
which,  however,  would  make  the  scheme  expensive  to  apply.  This  would  also 
means  waste  of  computational  effort  since  it  would  not  change  the  situation  of 
well  resolved  structures  at  the  smooth  regions  of  the  solution.  Moreover,  the 
hxed  order  of  the  hnite  diherences  does  not  provide  the  exponential  resolution 
of  spectral  schemes. 

This  article  aims  at  the  conjugation  of  the  spectral  and  WENO  methods  when 
solving  hyperbolic  equations  with  discontinuous  solutions.  The  general  idea  is 
to  build  a  multi-domain  scheme,  forming  a  global  adaptive  mesh  composed 
of  WENO  and  spectral  subdomains  in  a  way  that  discontinuities  of  the  solu¬ 
tion  are  always  contained  within  WENO  subdomains  and  the  smooth  com¬ 
ponents  remain  in  spectral  ones.  A  smoothness  measurement  device  triggers 
the  switching  of  the  subdomains  from  spectral  to  WENO  and  reciprocally,  ac¬ 
cording  to  the  local  behavior  of  the  solution.  For  instance,  the  algorithm  takes 
care  of  moving  discontinuities  by  changing  the  spectral  subdomains  in  their 
way  to  WENO  ones  and  switching  WENO  subdomains  in  their  trails  to  spec¬ 
tral  subdomains,  guaranteeing  higher  numerical  efficiency  than  the  classical 
single-domain  WENO  method  and  non-oscillatory  solutions  when  compared 
to  the  classical  single-domain  spectral  method.  Moreover,  the  application  of 
spectral  methods  at  smooth  regions  avoids  the  heavy  machinery  employed  by 
the  characteristic-wise  WENO  hnite  difference  algorithm  such  as  the  evalua¬ 
tion  of  the  Jacobian  of  the  huxes,  the  global  Lax-Friedrichs  hux  splitting  and 
the  forward  and  backward  characteristics  projections. 

The  multi-domain  hybrid  Spectral- WENO  method  (Hybrid)  here  proposed 
can  also  be  thought  of  as  an  improvement  to  the  classical  spectral  method 
by  using  the  ’’WENO  technique”  to  approximate  discontinuities,  in  the  same 
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way  that  physical  hltering  is  used.  Based  on  the  multi-domain  framework, 
the  Gibbs  phenomenon  is  effectively  avoided,  since  there  will  be  no  spectral 
approximation  of  discontinuities,  which  also  discards  the  need  of  any  postpro¬ 
cessing  technique. 

This  paper  is  organized  as  follows:  Approximation  theory  of  Spectral  meth¬ 
ods  is  briefly  discussed  in  Section  2.  In  Section  3,  WENO  Finite  Differences 
schemes  are  described  and  Harten’s  multi-resolution  algorithm  is  presented  in 
Section  4.  The  Hybrid  method  along  with  its  interfaces  treatment  and  subdo¬ 
main  switching  procedure  are  introduced  in  Section  5.  In  Section  6,  we  apply 
the  Hybrid  method  to  the  standard  shock-tube  problems  with  Riemann  ini¬ 
tial  data,  to  the  Shock-Entropy  wave  interaction  and  the  Blastwave  problems. 
Conclusions  are  given  in  Section  7. 


2  Spectral  Methods 


We  start  this  section  with  a  quick  description  of  Galerkin  and  collocation  spec¬ 
tral  methods,  which  in  their  simplest  forms,  make  use  of  a  global  smooth  basis 
{0fc(a;),  /c  =  0, . . . ,  N}  to  represent  the  numerical  solution.  After  that,  we  dis¬ 
cuss  the  Gibbs  Phenomenon  for  piecewise  smooth  functions  and  the  utilization 
of  hltering  to  stabilize  the  spectral  discretization  of  hyperbolic  conservation 
laws. 


2.1  Spectral  Collocation  Methods 


In  the  Galerkin  approach,  the  function  f{x)  is  projected  into  the  set  of  basis 
function  0fc(a;)  as 

N 

PNf{x)  =  '^ak(l)k{x),  (1) 

k=0 

where  is  the  projection  operator.  The  coefficients  are 


f{x)(l)k{x)w{x)dx, 


(2) 


with  the  appropriate  weight  function  w{x)  that  depends  on  the  nature  of  the 
physical  solution: 

•  For  periodic  problems,  the  natural  basis  functions  are  the  trigonometric 
polynomials  of  degree  k,(j)k{x)  =  with  weight  function  w{x)  =  1  and 
a;  e  [-1,1). 
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•  For  non-periodic  problems  in  a  finite  domain  x  G  [—1, 1],  Chebyshev  poly- 

_  1 

nomials  (pk{x)  =  Tk{x)  with  w{x)  =  (1  —  x'^)  2  or,  Legendre  polynomials 
4>k{x)  =  Lk{x)  with  w{x)  =  1  are  the  most  used  bases. 

The  collocation  version  of  spectral  methods  approximates  a  function  f{x)  by 
an  interpolating  polynomial  given  by 

N  N 

Infix)  =  ttkfkix),  ttk  =  Y^ifi^i)fk{Xi),  (3) 

fc=0  i=0 

where  is  the  interpolation  operator,  Xi  and  cu*  are  the  Gauss-Lobatto 
quadrature  nodes  and  weights  respectively  (Gauss-Radau  and  Gauss  nodes 
can  also  be  used).  Alternatively, 


N 

Infix)  =  Yfi^j)9jix), 

j=0 


(4) 


where  Qjix)  are  the  cardinal  functions:  (trigonometric)  Lagrangian  interpo¬ 
lation  polynomials  of  degree  N  such  that  Qjixi)  =  Sij.  For  the  Lagrangian 
interpolation  polynomials  based  on  the  Ghebyshev  Gauss-Lobatto  points  Xi  = 
cosini /N),i  =  0, . . . ,  N,  the  cardinal  functions  are  given  by 


CjN‘^ix  —  Xj) 


(5) 


where  Cq  =  Cjv  =  2,  =  l,j  =  1, . . . ,  N  —  1  and  Tj^ix)  is  the  A^_th  de¬ 

gree  Ghebyshev  polynomial  of  the  first  kind.  The  derivatives  of  fix)  at  the 
collocation  points  Xi  can  be  computed  via  equations  (3)  or  (4).  The  former 
makes  use  of  the  fast  cosine  transform  (GFT)  algorithm  and  the  latter  uses  a 
matrix- vector  algorithm.  More  details  can  be  found  in  [1]. 


2.2  Conservation  Laws:  Gibbs  Phenomenon  and  Filtering 


The  approximation  error  in  spectral  methods  depends  only  on  the  regularity 
of  the  approximated  function.  A  typical  error  estimate  is  of  the  form 

1 

\f{x)  -  P„f{x)\  <  cm-r  {^£'  l/^Kjpdf)  "  ,  (6) 

where  G  is  a  constant  independent  of  N  and  fl’'^  denotes  the  p_th  derivative  of 
/.  We  see  that  the  approximation  error  decays  as  0(A^“^)  for  any  function 
and  if  the  function  is  analytic  then 

|/(a;)-P,/(a;)|  <Ge-“^,  (7) 
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for  some  a  >  0,  resulting  in  exponential  convergence.  Similar  results  hold  for 
the  pseudospectral  (collocation)  formulation.  However,  in  the  case  of  piecewise 
smooth  functions,  the  order  of  accuracy  is  reduced  to  0(1)  due  to  the  well 
known  Gibbs  phenomenon. 

Figure  2  shows  the  Fourier  collocation  approximations  to  a  sawtooth  function 
with  a  discontinuity  at  a;  =  0.  Note  that  the  over-  and  under-shoot  oscillations 
do  not  decrease  with  the  increasing  number  of  grid  points  N. 


Fig.  2.  Gibbs  phenomenon  in  the  spectral  approximation  of  a  sawtooth  function. 


Other  important  features  of  spectral  methods  are  their  inherent  non-dissipativity 
and  non- dispersiveness,  which  are  good  properties  when  important  conserva¬ 
tion  principles  need  to  be  respected  at  the  numerical  level.  This  lack  of  dissi¬ 
pation,  however,  becomes  an  issue  when  dealing  with  discontinuous  solutions 
of  conservation  laws,  in  particular,  when  spectral  methods  are  applied  to  non¬ 
linear  hyperbolic  equations  in  the  conservation  form  and  the  problem  of  an 
entropy  satisfying  solution  arises.  In  fact,  there  is  no  artihcial  dissipation  in 
spectral  methods  to  indicate  that  their  solutions  are  limits  of  a  dissipative 
process.  This  problem  had  been  addressed  in  [12]  and  the  references  contained 
therein,  where  it  has  been  shown  that  with  a  suitable  addition  of  a  spectrally 
small  dissipation  to  the  high  modes,  the  method  is  stable  and  entropic  solu¬ 
tions  are  obtained. 

The  utilization  of  low  pass  hltering  in  spectral  methods  helps  both  to  mitigate 
the  Gibbs  phenomenon  and  to  introduce  the  necessary  artihcial  dissipation  for 
hyperbolic  problems.  Gonsidering  the  Fourier  approximation  (The  process  is 
analogous  for  the  Ghebyshev  case): 

N 

U{x)=  a;e[-l,l),  (8) 

k=-N 

hltering  is  introduced  by  means  of  the  modihed  sum 

/;w=  E  (-I) (9) 

k=-N  / 
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where,  following  Vandeven  [13],  cr(a;)  is  a  p_th  order,  p  >  1,C^[— 1,1],  filter 
function  satisfying 


a(0)  =  l,  a(±l)  =  0, 
ad)(0)  =  0,  a(^)(±l)  =  0,  j  <  p, 


(10) 


where  denotes  its  j_th  derivative. 

It  has  been  shown  in  [13]  that  the  hltered  sum  (9)  approximates  a  discontin¬ 
uous  function  in  smooth  regions  with  exponential  accuracy. 

In  this  work  we  make  use  of  the  Exponential  hlter  function 

cr(a;)  =  exp  a|a;|^^  ,  (11) 

where  —1  <  uj  =  k/N  <  1,  |/c|  =  0, . . . ,  iV  ,/?  is  the  order  of  the  hlter  and 
a  =  —  Ine,  where  e  is  the  machine  zero.  In  hgure  3  the  solution  of  the  period¬ 
ical  inviscid  Burgers  equation  is  obtained  with  a  Fourier  collocation  method 
stabilized  with  a  16_th  order  Exponential  hlter. 


Fig.  3.  (Left)Exact  and  Fourier  spectral  solution  of  the  Burgers  equation  with  a 
16_th  order  Exponential  filter  and  N  =  128  collocation  points.  (Right)  Convergence 
of  the  solution  with  N  =  128,  256, 512.  The  error  is  in  the  loglO  scale. 

It  is  worth  saying  that  no  amount  of  hitering  other  than  a  heavy  one  (p  <  2), 
can  remove  the  Gibbs  oscillations  from  the  solution.  The  use  of  such  a  hlter, 
however,  also  incurs  in  removal  of  small  and  median  scale  structures  rendering 
the  solution  unusable  and  inaccurate. 

In  practice,  when  solving  diherential  equations  one  uses  a  high  order  Expo¬ 
nential  hlter  at  every  time  step  to  maintain  stability  and  a  post-processing 
technique  [4,27,29]  at  the  end  of  the  calculation  to  recover  a  non-oscillatory 
solution  from  the  oscillatory  one.  These  post-processing  techniques,  however, 
are  not  practical,  particularly,  for  multi-dimensional  problems.  As  it  is  evident 
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from  figure  3,  in  order  to  recover  the  full  accuracy  in  any  region  where  the 
function  is  continuous,  one  has  to  use  a  different  idea.  In  the  next  section  we 
introduce  non-oscillatory  hnite  differences  schemes  as  the  next  step  to  build 
this  new  idea  in  the  form  of  a  hybrid  scheme. 


3  Weighted  essentially  non-oscillatory  schemes 


In  this  section  we  describe  essentially  non-oscillatory  conservative  hnite  differ¬ 
ence  schemes  to  be  applied  at  systems  of  conservation  laws  described  by  the 
general  equation: 

ut  +  f{u),  =  0.  (12) 

We  will  hrst  present  a  general  high  order  conservative  central  scheme  approx¬ 
imation  to  (12).  After  that,  we  will  describe  the  adaptive  stencil  choosing 
process  of  ENO  and  WENO  schemes  that  leads  to  their  non-oscillatory  prop¬ 
erties. 


3.1  Conservative  Schemes 


Being  conservative  is  an  important  property  of  a  numerical  approximation  due 
to  the  Lax-Wendroff  theorem  which  states  that  numerical  solutions  of  conser¬ 
vative  schemes,  whenever  they  converge,  they  converge  to  weak  solntions  of 
the  conservation  law  [9]. A  conservative  hnite-difference  formulation  for  hyper¬ 
bolic  conservation  laws  reqnires  high-order  consistent  nnmerical  huxes  at  the 
cell  bonndaries  in  order  to  form  the  flux  difference  across  the  nniformly-spaced 
cells. 


Consider  an  uniform  grid  dehned  by  the  points  Xi  =  iAx,  i  =  0, . . . ,  N,  which 

are  also  called  cell  centers,  with  cell  boundaries  given  by  x.  i  =  Xi  ± 

1+2 

Eqnation  (12)  is  semi-discretized  in  time  by  the  method  of  lines,  yielding  the 
following  system  of  ordinary  differential  equations 


duiit) 

dt 


dx 


=  0,...,A^, 


X=Xi 


(13) 


where  Ui{t)  is  a  numerical  approximation  to  the  cell  center  valne  u{xi,t)- 


The  conservative  property  of  the  spatial  discretization  is  obtained  by  implicitly 
dehning  a  numerical  flux  function  h{x)  as 


fix) 


1 


(14) 


such  that  the  spatial  derivative  in  (13)  is  exactly  approximated  by  a  conser¬ 
vative  hnite  difference  formula  at  the  cell  boundaries, 


duiit) 

dt 


(15) 


where  h.  i  =  h{x  ±  ^).  High  order  polynomial  interpolations  to  h.  i  are 
i±  2  2 

computed  using  known  grid  values  of  /. 


Note  that  the  order  of  the  polynomial  interpolation  of  h.  i  will  set  the  order 

i±  2 

of  the  spatial  approximation  of  the  overall  scheme  (15).  Due  to  the  symmetry 
of  the  stencils  used  for  the  polynomial  interpolations,  the  order  of  approxima¬ 
tion  of  the  right-hand  side  of  (15)  will  be  the  same  order  of  the  polynomial 
approximation  of  h{x),  even  after  division  by  Ax,  and  not,  as  expected,  one 
order  less. 


Let  us  now  solve  the  problem  of  hnding  a  kAh.  order  central  polynomial  ap¬ 
proximation  to  h.  1.  Let 
*+2 


h{x)  =  Oo  -l-  aix  -1-  ■  ■  ■  -|-  ttk-ix^  ^  (16) 

be  a  (A;  —  l)_th  degree  polynomial  approximation  to  h{x)  and  f{x)  be  the 
polynomial  (also  of  degree  k  —  1)  obtained  by  integration  of  (14)  after  the 
substitution  of  h{x)  by  h{x).  The  polynomial  coefficients  are  found  after  eval¬ 
uating  f{x)  at  the  grid  points  of  any  k  nodes  stencil  around  For  in¬ 

stance,  if  A;  =  3,  the  coefficients  oq,  01,02  can  be  obtained  as  functions  of  the 
grid  point  values  and  /i+i.  Once  the  coefficients  are  fonnd,  the  valnes 

/  1  =  h  i  +  0{Ax^)  can  be  compnted  and  snbstituted  into  (15). 

2  ^^2 

The  scheme  above  is  centralized  at  a;  1  only  if  k  is  even,  otherwise,  it  is  shifted 

*+2 

a  half-cell  to  the  left  (npstream)  or  to  the  right  (downstream).  One  often 
takes  the  upstream  version  of  an  odd  order  central  scheme  dne  to  the  inherent 
dissipation  of  upstream  schemes,  necessary  to  shock-capturing.  For  the  sake 
of  completeness,  let  us  now  take  the  case  k  =  5  and  make  the  computations 
to  obtain  the  5_th  order  approximation  to  the  spatial  derivative  in  (15). 

Thus,  after  substitution  of  a  4_th  order  polynomial  in  (14)  and  integration, 
we  obtain 


f{x)  =  ao+aiX+a2  x^  + 


Ax^ 

~i2 


I  +0,3  X  + 


xAx'^ 


I  +(^4  X  + 


x‘^Ax‘^ 


,  (17) 


which,  if  compnted  at  the  nodes  {xi-2,  ■  ■  ■ ,  XiJ^2}  will  determine  the  coefficients 
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(18) 


{ao, . . . ,  04}  in  terms  of  {/i_2,  •  •  • ,  fi+2},  yielding 


^  (2/i-2  —  13/i-i  +  4:7 fi  +  27/j+i  —  3/4+2) , 
*+2  oU 

and,  analogously, 

/,  1=4  (2/*-3  -  13/4-2  +  47/4-1  +  27/4  -  3/4+1) . 
*“2  dU 


Both  are  5_th  order  approximations  to  h.  1  and  h.  1,  respectively: 

*+  2  *”2 


d^f 
60  dx^ 


+  0(Aa;®). 

X=Xi 


(19) 


(20) 


The  leading  error  order  is  the  same  for  both  stencils,  so  they  will  cancel  out 
after  substitution  of  (18)  and  (19)  at  (15),  yielding  a  5_th  order  approximation 
as  pointed  out  beforehand. 


3.2  WENO  Convex  Combination  of  Stencils 


The  central  hnite  difference  scheme  presented  above  will  suffer  from  spurious 
oscillations  if  discontinuities  are  inside  the  interpolating  stencils.  The  main 
idea  of  ENO  schemes  is  to  bias  as  much  as  necessary  the  central  scheme  to 
avoid  the  inclusion  of  the  shock  inside  the  stencil.  For  instance,  a  k-th  order 
ENO  approximation  assigns  smoothness  weights  to  all  k  nodes  stencils  around 
the  interpolating  point  and  chooses  the  smoothest  one  to  perform  k-th  order 
polynomial  interpolations  of  the  numerical  flux. 

Let  us  identify  a  particular  stencil  by  its  left-shift  r.  We  have  a  polynomial 

interpolation  to  h,  1  for  each  r,  given  by 
*+2 

/r(a^4+i )  =  14  Crjfi-r+j  =  h{x  l)  +  0{Ax’^),  (21) 

2  j=o  2 

where  the  Crj  are  the  Lagrangian  interpolation  coefficients  ([24]).  The  Crj  de¬ 
pend  on  the  order  k  of  the  approximation  and  also  on  the  left-shift  parameter  r, 
but  not  on  the  values  fi.  Since  r  can  vary  from  0  to  A;  —  1,  we  have  k  distinct 
interpolating  polynomials  to  choose  from;  all  of  them  yielding  a  /cTh  order 
approximation,  once  /  is  smooth  inside  the  intervals  considered.  The  above 
process  is  called  the  reconstruction  step,  for  it  reconstructs  the  values  of  h{x) 

at  the  cell  boundaries  of  the  interval  Ii  =  \x.  i,x.  i  ]  from  the  cell  average  val- 

i  2  2 

ues  f{x)  in  the  intervals  Sr  =  {Uj=o  h-r+jC  =  0, . . . ,  /c  —  1}  =  {xi^r,  •••, 
with  s  =  k  —  r  —  1. 
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However,  at  smooth  regions,  the  collection  of  all  k  stencils  carry  information 

for  an  approximation  of  order  higher  than  k.  The  Weighted  Essentially  Non- 

Oscillatory  scheme  (WENO)  is  an  improvement  over  ENO  for  it  uses  a  convex 

combination  of  all  available  polynomials  for  a  hxed  k,  assigning  zero  weights 

to  stencils  containing  discontinuities.  This  yields  a  {2k  —  1)  order  scheme  at 

smooth  parts  of  the  solution.  The  general  WENO  flux  /,  i  is  dehned  as 

1+2 


k-1 


where 


with 


Cr 

(ir  +  /Sr)»' 


(23) 


Here,  Cr  are  the  ideal  weights  for  the  convex  combination,  the  ones  that  at 
the  absence  of  discontinuities  provide  a  (2r  —  l)_th  order  of  approximation  at 
(22)  ([24]);  e  =  10“^°  is  a  small  parameter  to  avoid  division  by  zero,  p  =  2  is 
chosen  to  increase  the  difference  of  scales  of  distinct  weights  at  non-smooth 
parts  of  the  solution  and  is  a  measure  of  the  smoothness  of  polynomials 
on  the  r_th  stencil; 


ISr 


dx. 


(24) 


When  the  interpolating  polynomial  on  a  given  stencil  is  smooth,  the  smooth¬ 
ness  indicator  ISr  is  relatively  much  smaller  than  those  of  stencils  where  the 
polynomial  has  discontinuities  in  its  first  k  —  1  derivatives.  Therefore,  discon¬ 
tinuous  stencils  receive  a  close  to  zero  weight,  ~  0,  and  a  non-oscillatory 
property  is  achieved. 


Remark  1  There  are  different  choiees  for  the  smoothness  indieator  I Sk  and 
Taylor  analysis  reveals  relevant  properties  of  eaeh  of  them.  The  first  one  was 
proposed  in  [24]  and  was  based  on  divided  differenees.  In  [6],  an  improvement 
over  (24)  was  proposed  in  order  to  inerease  the  accuracy  of  the  scheme  at 
critical  points  (zero  derivatives)  of  the  solution. 

For  a  system  of  Conservation  Laws  such  as  the  Euler  equations,  the  eigen¬ 
vectors  and  eigenvalues  of  the  Jacobian  of  the  flux  are  computed  via  the  Roe 
Average  method.  Global  Lax-Friedrichs  flux  splitting  is  used  to  split  the  flux 
into  its  positive  and  negative  components.  Artificial  dissipation  based  on  the 
modulus  of  the  eigenvalues  is  added  in  order  to  obtain  a  smoother  flux.  The 
resulting  positive  and  negative  flux  components  are  then  projected  into  the 
characteristic  helds  using  the  left  eigenvectors  to  form  the  positive  and  nega¬ 
tive  characteristic  variables  at  each  grid  cell  center.  Then,  high-order  WENO 
polynomial  reconstruction,  as  described  above,  is  used  to  obtain  these  char¬ 
acteristic  variables  components  at  the  grid  cell  boundaries,  which,  after  sum- 
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mation,  are  finally  projected  back  to  physical  space  via  the  right  eigenvectors. 
The  details  of  this  algorithm  can  be  found  in  [24],  These  characteristic  vari¬ 
ables  projections  are  the  expensive  part  of  the  WENO  scheme  when  applied 
to  systems  of  equations.  They  are  necessary  because  high  order  approxima¬ 
tion  is  not  achieved  within  the  framework  of  the  conservative  variables.  A 
hybrid  Spectral- WENO  scheme  would  save  substantial  computational  effort 
by  avoiding  the  characteristic  projections  at  all  smooth  parts  of  the  solution. 


4  Multi-Resolution  Analysis 


The  successful  implementation  of  the  Hybrid  method  depends  on  the  ability 
to  obtain  accurate  information  on  the  smoothness  of  a  function.  In  this  work, 
we  employ  the  Multi-Resolution  (MR)  algorithms  by  Harten  [22,23]  to  detect 
the  smooth  and  rough  parts  of  the  numerical  solution.  The  general  idea  is 
to  generate  a  coarser  grid  of  averages  of  the  point  values  of  a  function  and 
measure  the  differences  (MR  coefficients)  di  between  the  interpolated  values 
from  this  sub-grid  and  the  point  values  themselves.  A  tolerance  parameter  Smr 
is  chosen  in  order  to  classify  as  smooth  those  parts  of  the  function  that  can 
be  well  interpolated  by  the  averaged  function  and  as  rough  those  where  the 
differences  di  are  larger  than  the  parameter  eMR-  We  shall  see  that  the  order  of 
interpolation  is  relevant  and  the  ratio  between  di  of  distinct  orders  may  also 
be  taken  as  an  indication  of  smoothness. 


Let  us  start  by  showing  two  examples  where  one  can  notice  the  detection 
capabilities  of  the  Multi-Resolution  analysis  that  will  be  presented  below. 
The  left  and  right  hgures  of  hgure  4  show  the  piecewise  analytic  function 


f{x)  = 


10  -I-  — 1  <  a;  <  —0.5 


X 


-0.5  <  a;  <  0 


sin(27ra;)  0  <  a;  <  1 


(25) 


and  the  density  (p)  of  the  Mach  3  Shock-Entropy  wave  interaction  problem 
[24],  as  computed  by  the  classical  hfth  order  WENO  hnite  difference  scheme, 
respectively. 


The  test  function  (25)  has  a  jump  discontinuity  at  a;  =  —0.5  and  a  disconti¬ 
nuity  at  its  first  derivative  at  a;  =  0.  One  can  see  that  at  each  grid  point  the 
differences  di  decay  exponentially  to  zero  inside  the  analytical  pieces  of  the 
function  when  the  order  of  interpolation  increases  from  Umr  =  3  to  Umr  =  8. 
At  the  discontinuity  x  =  0.5,  the  measured  differences  di  are  0(1)  and  remain 
unchanged  despite  the  increase  of  the  interpolation  order.  Similar  behavior  is 
exhibited  at  the  derivative  discontinuity  at  a;  =  0  with  a  smaller  amplitude. 
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Fig.  4.  (Left)  The  third  and  eighth  order  MR  coefficients  di  of  the  piecewise  analytic 
function.  (Right)  The  third,  fifth  and  seventh  order  MR  coefficients  di  of  the  density 
f{x)  =  p  of  the  Mach  3  Shock-Entropy  wave  interaction  problem. 

Also,  in  the  right  hgure  of  hgure  4,  the  density  of  the  Mach  3  Shock-Entropy 
wave  interaction  problem  and  the  corresponding  MR  coefficients  di  are  shown 
for  the  third,  hfth  and  seventh  order  Multi-Resolution  analysis.  The  location 
of  the  main  shock  is  at  x  ~  2.73  and  the  shocklets  behind  the  main  shock  are 
well  captured.  The  high  frequencies  behind  the  main  shock  are  much  better 
distinguished  with  the  higher  orders. 

Averaging  a  function  corresponds  to  hlter  the  upper  half  of  the  spectrum.  The 
main  idea  of  Hartens  smoothness  classihcation  is  to  measure  how  distant  the 
actual  values  of  the  function  are  from  being  predicted  through  interpolation 
of  the  lower  half  of  the  frequencies  contained  in  the  sub-grid  of  averages. 
We  now  describe  a  detailed  construction  of  the  sub-grid  of  averages  and  its 
corresponding  interpolating  polynomial,  hnishing  with  a  worked  example. 

Given  an  initial  number  of  grid  points  Nq  and  grid  spacing  Axq,  consider  the 
set  of  nested  dyadic  grids  {G^,  0  <  k  <  L},  dehned  as: 


=  =  (26) 

where  xf  =  iAxk,  Axk  =  2^Axo,  Nk  =  2~^No.  For  each  level  A:  >  0  we  dehne 
the  set  of  cell  averages  {/f ,  i  =  1, . . . ,  Nk}  at  xf  of  a  function  /(x): 


fl 


1 

Axfc 


f{x)dx, 


(27) 


and  /°  =  /°.  Let  f^i-i  be  the  approximation  to  f^i-i  by  the  unique  polynomial 
of  degree  2s  that  interpolates  |/|  <  s  at  where  r  =  2s-|-  1  is  the  order 
of  approximation. 

The  approximation  differences,  also  called  multiresolution  coefficients,  d’^  = 
/2A1  —  f 21^11  the  k_th  grid  level  and  grid  point  x*,  have  the  property  that 
if  /(x)  has  p  —  1  continuous  derivatives  and  a  jump  discontinuity  at  its  p_th 
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derivative,  then 


4  «  (  ^  ^  (28) 
Axlfi  for  p>  r 

where  [■]  denotes  the  magnitude  of  the  jump  of  the  function  inside. 

From  formula  (28)  it  follows  that 

I4“^l  ~  2-P|df|,  where  p  =  min{p,r}.  (29) 

Equation  (29)  shows  that  away  from  discontinuities,  the  MR  coefficients  df 
diminish  in  size  with  the  rehnement  of  the  grid;  close  to  discontinuities,  they 
remain  the  same  size,  independent  of  k.  The  MR  coefficients  df  were  used  in 
[23]  in  two  ways.  First,  hner  grid  data  /?  were  mapped  to  its  M  level  mul¬ 
tiresolution  representation  =  (d),  •  •  •  ,df,  //“)  to  form  a  multiscale  version 
of  a  particular  scheme,  where  truncation  of  small  quantities  with  respect  to  a 
tolerance  parameter  decreased  the  number  of  flux  computations.  Secondly,  the 
MR  coefficients  df  also  acted  as  a  shock  detection  mechanism  and  an  adaptive 
method  was  designed  where  a  second-order  Lax-Wendroff  scheme  was  locally 
switched  to  a  hrst-order  accurate  TVD  Roe  scheme,  whenever  dj  was  bigger 
than  ^ M R' 

Equation  (28)  also  indicates  that  the  variation  of  the  MR  order,  Umr,  can  give 
additional  information  on  the  type  of  the  discontinuity.  Nevertheless,  in  this 
work,  we  will  be  limited  at  using  only  the  hrst  level  k  =  lof  the  multiresolution 
coefficients  and  we  shall  drop  the  superscript  1  from  the  d}  from  here  on  unless 
noted  otherwise. 

Hence,  to  hnd  di,  the  idea  is  to  construct  a  piecewise  polynomial  Pk{x)  of 
degree  k  =  Umr  using  k  +  1  computed  average  values  of  /j,  /j,  at  the  equi- 
spaced  grid  Xi  such  that 


Pkix^)  =  fix,)  +  OiAx^^^),  (30) 

and 

di  =  fi  -  Pk{xi).  (31) 

Given  a  tolerance  level  eMR,  the  smoothness  of  the  function  f{x)  at  Xi  would 
then  be  checked  against  the  magnitude  of  the  di,  namely: 

{Mil  <  Smr  ^  solution  is  smooth. 

(32) 

Mil  >  ^mr  ^  solution  is  non-smooth. 

The  algorithm  for  computing  the  MR  coefficients  di  is  given  next. 
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4-1  Computing  the  MR  Coefficients 


Consider  an  equi-spaced  grid  {xi  =  iAx,  i  =  —m, . . .  ,0, . . . ,  N, . . . ,  N  +  M} 
where  Ax  is  the  constant  grid  spacing.  N  can  be  an  odd  or  an  even  number. 
Depending  on  N  and  the  even  or  odd  order  of  the  MR  Analysis  Umr,  the 
number  of  ghost  points  m  and  M  required  are  given  in  table  I. 


N 

Umr 

m 

M 

odd 

odd 

Umr  +  1 

Umr  +  1 

even 

odd 

Umr  +  1 

nMR 

odd 

even 

Umr 

Umr  +  2 

even 

even 

Umr 

Umr  +  1 

Table  I 

The  number  of  ghost  points  m  and  M  required  for  the  MR  Analysis. 


Given  the  grid  point  values  of  the  function  f{x),  the  average  values  are  com¬ 
puted  as 


h 


{f2i  +  f2i+l)  , 


m  N  +  M  —  1 
¥’■■■’  2 


(33) 


We  construct  a  piecewise  k  =  Umr  degree  polynomial  Pk{x)  using  the  k  +  1 
computed  average  values  of  the  given  function,  /j  such  that 


Pkixi)  =  f{xi)  +  0(Aa;^+^). 


(34) 


The  polynomial  Pk{xi),  I  =  \m  and  L  =  /  —  1  or  L  =  /  if  /c  is  odd  or  even, 
respectively,  can  be  written  as 


i+L 

Pk{Xi)  =  arfr-  (35) 

r=i—l 

However,  since  the  coefficients  a  depend  only  on  Xi  and  do  not  depend  on  the 
function  f{x),  the  Pk{xi)  can  be  written  as: 


Pk  i,Xi) 


Yir=-l  Oirfi+r,  mod(b  2)  =  0 
Y.r=-iPr+ifi+T,  mod(b2)  =  1 


(36) 


In  the  case  of  mod(q  2)  =  1, 


jS-r  =  (Xr,  r  =  —l,...,L.  (37) 

Furthermore,  if  Umr  is  even,  the  coefficients  a  are  symmetric  about  r  =  0, 
namely,  a-r  =  Or,  r  =  1, . . . ,  L. 
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The  desired  coefficients  a  are  computed  by  requiring  Pk{x)  to  be  equal  to 
each  of  the  first  k  +  1  monomials  f{x)  =  l,x,x‘^, . . . ,  x^  and  evaluated  at  any 
grid  point  x  =  x*.  For  simplicity,  we  take  x*  =  0.  The  fi  are  evaluated  for 
i  =  —I, . . . ,  L.  This  procedure  results  in  a  system  of  linear  equations,  Aa  =  6, 
where 


A  = 


l'  1  ...  1 

-2/  +  (-2/  +  1)  ...  2L  +  (2L  +  1) 


\  (-2/)^  +  {-21  +  1)^  . . .  (2L)^  +  (2L  +  1)^  J 

and  A  is  a  matrix  of  size  (L  +  /  +  1)  x  (L  +  /  +  1) 


l\ 

Oi-l+l 

0 

,  "cr  = 

,  b  = 

\  J 

vO; 

(38) 


Using  (36),  the  k-th  order  Multi- Resolution  coefficients  di  at  Xi  can  be  com¬ 
puted  as 


di  =  fi-Pk{xi)  f  =  0,...,A. 


(39) 


One  can  also  evaluate  the  a  by  matching  the  terms  in  the  Taylor  series  ex¬ 
pansion  using  (34)  and  (35)  to  any  desired  order,  however  this  procedure  may 
become  cumbersome  for  high  order  k. 


Example 


To  illustrate  the  procedure  above,  we  will  construct  two  unique  local  poly¬ 
nomials  with  k  =  Umr  =  3,  such  that  Pk{xo)  =  f{xo)  +  0{Ax^~^^)  and 
Pk{xi)  =  f{x^)  +  0{Ax’^^^). 

To  construct  the  desired  polynomials  one  needs  to  hnd  the  unique  coefficients 
{q;_2,  a-i,  tto,  tti}  and  /do, /5i, /52}  such  that 

a_2/-2  +  a-i/-i  +  CKo/o  +  ai/i  =  f{xo)  +  0(Aa;'‘)  (40) 

and 

P-if-i  +  Pofo  +  Pifi  +  ^2/2  =  f{xi)  +  0{Ax^).  (41) 
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The  system  of  equations,  (38),  becomes 


(  1  1  ll\ 

^l\ 

-7-3  15 

CX—l 

0 

,  7?  = 

,  b  = 

25  5  1  13 

ao 

0 

^-91  -9  1  35  y 

[  a,  ) 

loj 

Solving  this  system  yields 


OL-2  —  — 


3 
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55 

“““S’ 


and  {/5_i  =  q;i,/3o  =  q;o,/9i  =  a-i,/92  =  tt-2}- 


(42) 


(43) 


Remark  2  The  tolerance  parameter  Cmr  determines  the  dynamic  activation 
of  the  spectral  and  WENO  spatial  discretizations  along  the  various  subdomains 
of  the  hybrid  method.  While  a  too  small  value  of  Cmr  activates  the  more  ex¬ 
pensive  WENO  method  at  subdomains  where  the  solution  is  smooth,  a  larger 
value  activates  the  spectral  method  at  a  subdomain  with  low  spatial  resolution, 
generating  oscillations.  Cmr  also  bears  a  straight  relation  with  the  interpola¬ 
tion  order  n MR-  High  Umr  values  decrease  the  size  of  Cmr  one  needs  to  chose, 
since  high  freguencies  are  less  mistaken  by  gradient  jumps.  The  general  guide¬ 
line  is  to  start  with  a  value  for  Umr  at  least  egual  to  the  order  of  the  WENO 
method  and  increase  it  according  to  the  complexity  of  the  solution.  Eor  in¬ 
stance,  Umr  =  5  is  a  good  choice  for  the  piecewise  smooth  solution  of  the  SOD 
problem,  the  Entropy  problem  would  work  better  with  Umr  =  7.  Eor  most  of 
the  flows  with  shock  that  were  tested,  the  value  of  Cmr  =  10“^  yielded  a  good 
balance  between  computational  speed  and  accuracy  of  the  numerical  solution. 


5  The  Multi-Domain  Hybrid  Spectral— WENO  Method 


The  main  idea  is  straightforward: 


Partition  the  physical  domain  into  a  number  of  equal  sized  subdomains,  avoid 
Gibbs  phenomenon  by  treating  discontinuities  with  shock  capturing  non-oscillatory 
WENO  methods  and  increase  the  numerical  efficiency  by  treating  the  smooth 
parts  of  the  solution  in  with  spectral  methods. 


In  the  case  of  a  stationary  discontinuity,  one  only  needs  appropriate  interface 
conditions  to  transmit  data  between  the  subdomains.  A  moving  shock  situ¬ 
ation  requires  in  addition  a  smoothness  measurement  routine  to  keep  track 
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of  the  high  gradients  but  also  a  subdomain  switching  algorithm  in  order  to 
change  spectral  domains  to  WENO  ones,  once  the  shock  gets  closer  to  them, 
and  switch  back  WENO  subdomains  to  spectral  ones,  once  the  shock  leaves 
them.  These  same  tools  will  do  the  work  in  the  case  when  discontinuities  are 
being  created  in  spectral  subdomains,  as  we  shall  see  later  in  Section  6,  when 
we  simulate  the  interaction  of  a  Mach  3  shock  with  a  sinusoidal  density  wave. 
In  this  section  we  describe  the  necessary  treatment  of  the  interfaces  between 
the  subdomains  and  the  algorithm  used  to  switch  the  type  of  a  subdomain, 
from  WENO  to  spectral  and  vice  versa. 


5.1  Interfaces  treatment 


One  of  the  main  components  of  the  Hybrid  scheme  is  the  exchange  of  infor¬ 
mation  at  the  interfaces  between  subdomains  that  might  be  of  distinct  types. 
The  correct  setting  of  these  interfaces  is  of  great  relevance  in  order  to  avoid 
loss  of  accuracy  or  numerical  oscillations  that  could  contaminate  the  solution. 
The  Hybrid  scheme  deals  with  three  types  of  such  interfaces,  namely, 

(1)  Spectral-Spectral  interface; 

(2)  Spectral-WENO  interface; 

(3)  WENO-WENO  interface. 

We  shall  denote  the  grid  points  of  the  left  and  right  subdomains  by  x  and  y, 
respectively.  The  functional  values  at  these  points  are  denoted  by  s  and  w  for 
spectral  and  WENO  subdomains,  respectively.  Ng  and  are  the  number  of 
Chebyshev  collocation  points  of  the  spectral  grid  and  the  number  of  uniformly 
spaced  grid  points  of  the  WENO  subdomain,  respectively. 


Spectral-Spectral  Interface 


At  a  Spectral-Spectral  interface  the  solution  is  assumed  to  be  smooth,  oth¬ 
erwise  the  switching  algorithm  presented  in  the  next  subsection  would  have 
changed  the  subdomains  to  a  WENO  one.  This  fact,  along  with  the  high  order 
accuracy  of  the  spectral  solutions  in  each  subdomain  k,  ensures  that  a  simple 
average  of  the  two  functional  values  is  sufficient  for  obtaining  continuity  of  the 
solution  across  the  interface  for  the  Euler  Equations,  that  is, 

So  =  =  i  (so  + 

Other  types  of  interface  conditions,  such  as  the  imposition  of  the  solution  of 
the  Riemann  problem  at  the  spectral-spectral  interface,  yielded  no  discernible 
differences  in  the  cases  tested. 
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Spectral-WENO  interface 


At  a  Spectral-WENO  (WENO-Spectral)  interface,  the  ghost  points  of  the 
WENO  grid  {y-r,  ■  ■  ■ ,  V-i}-,  where  r  is  the  number  of  WENO  ghost  points,  in 
the  subdomain  k,  are  inside  the  spectral  subdomain  k  —  1;  and  the  end  point 
of  the  spectral  subdomain  k  —  1,  xq,  is  in  between  the  hrst  ghost  point  |/_i 
and  the  hrst  interior  point  yo  of  the  WENO  grid  in  subdomain  k.  The  spectral 
point  and  all  the  WENO  ghost  points  need  to  be  adjusted  for  continuity  of 
the  solution  at  the  interface.  The  idea  is  to  use  information  of  the  spectral 
subdomain  to  interpolate  at  the  ghost  points  and,  after  that,  to  interpolate 
the  spectral  boundary  point  using  the  surrounding  WENO  grid  points,  as 
detailed  here: 

•  The  solution  values  at  the  ghost  points  of  the  WENO  subdomain  are  com¬ 
puted  via  spectral  interpolation  using  the  data  (s(a;o), . . . ,  s(a;jv^))  of  the 
spectral  subdomain.  That  is, 

Ng 

^iVi)  =  Y.  i  =  -r, . . . ,  -1,  (45) 

j=0 

where  gj{x)  is  the  Lagrangian  interpolation  polynomial  of  degree  Ng  (see 
equation  5). 

•  The  functional  value  of  the  spectral  grid  at  xo  is  obtained  through  a  poly¬ 
nomial  interpolation  using  the  functional  values  of  the  interior  points  and 
the  updated  ghost  points  of  the  WENO  subdomain.  The  degree  of  the  in¬ 
terpolating  polynomial  should  match  the  order  of  the  WENO  method. 

The  hierarchy  of  interpolation  above  does  matter.  We  first  generate  the  ghost 
points  of  the  WENO  grid  because  these  are  inside  the  spectral  subdomain 
and,  therefore,  must  agree  with  the  spectral  solution.  Generating  the  spectral 
endpoint  Xq  at  first  would  use  wrong  information  from  the  ghost  points.  The 
case  of  a  WENO-Spectral  interface  is  completely  analogous  and  one  should 
also  respect  the  same  interpolation  hierarchy  as  above:  ghost  points  first. 


WENO-WENO  Interface 


It  will  be  through  the  interface  between  WENO  subdomains  where  possible 
discontinuities  in  the  solution  are  going  to  be  transmitted.  In  [16],  numerical 
experiments  show  that  WENO  or  Lagrangian  interpolation  are  not  conserva¬ 
tive  if  the  neighboring  subdomains  have  different  grid  spacing  Ax  and  that, 
once  this  condition  is  satisfied,  both  interpolations  are  conservative  and  be¬ 
have  equivalently.  In  this  work,  we  only  consider  WENO  subdomains  with  the 
same  grid  spacing  Ax  in  order  to  avoid  such  complications.  This  restriction 
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also  makes  ghost  points  of  a  subdomain  to  coincide  with  interior  points  of  the 
neighboring  subdomain  and,  for  practical  matters,  both  subdomains  can  be 
treated  as  a  single  larger  subdomain  with  no  interface,  allowing  shocks  to  pass 
through  as  if  they  were  ’’walking”  at  interior  points  of  a  WENO  discretization. 
If  we  denote  by  and  the  functional  values  of  the  left  and  right  WENO 
subdomains,  respectively,  we  have: 


^N^+i 

= 

=  ^Nw-i 

,i  =  1, 


,r, 


(46) 


where  and  the  size  of  the  adjacent  WENO  subdomains  are  assumed  to  be 
the  same. 


Example 


Let  us  now  show  by  means  of  a  numerical  example  that  the  hybrid  method  with 
the  above  interface  conditions  is  an  improvement  over  the  classical  WENO 
scheme.  Consider  the  inviscid  Burgers  equation 

I  u  (x,  0)  =  0.3  +  0.7  sin  {'kx)  , 

with  periodic  boundary  condition.  The  initial  conhguration  will  quickly  evolve 
into  a  moving  shock,  however,  at  this  point  we  will  analyze  the  pre-shock  error 
in  order  to  show  that  the  superior  accuracy  of  the  spectral  scheme  improves 
the  WENO  error  as  we  see  in  hgure  5.  Eor  the  classical  WENO  method  we 
discretize  the  interval  [—1,1]  with  400  points.  Eor  the  Hybrid  method,  the 
interval  is  subdivided  into  three  equal  sized  subdomains,  where  the  middle 
one,  [—0.3,  0.3],  is  a  spectral  subdomain.  The  subdomains  conhguration  is  kept 
hxed,  since  no  shock  will  arise  during  the  time  interval  of  the  experiment. 

The  number  of  points  in  the  WENO  subdomains  are  the  same  as  the  corre¬ 
sponding  regions  of  the  classical  WENO,  so  any  increase  of  accuracy  must  be 
due  to  the  spectral  subdomain.  Note  that  not  only  the  error  at  the  spectral 
subdomain  is  smaller  than  the  one  of  the  classical  WENO,  but  it  also  forces 
the  decreasing  of  the  overall  error,  showing  its  superior  performance  and  also 
the  conservation  properties  of  the  interface  conditions.  Another  aspect  that 
should  be  noticed  is  the  globality  of  the  spectral  scheme  that  can  be  observed 
when  looking  to  the  uniformity  of  the  error  at  the  spectral  subdomain. 
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Fig.  5.  (Dashed  line)  Inviscid  Burgers  Equation  solution;  (Continuous  line)  Error 
of  the  fifth-order  classical  WENO  with  400  points;  (Symbols)  Error  of  the  Hybrid 
method. 

5.2  The  Switching  Algorithm 


In  this  section  we  describe  the  snbdomain  switching  algorithm.  Since  we  are 
interested  in  solving  unsteady  moving  shocked  flows,  the  subdomains  must  be 
able  to  switch  from  one  type  of  subdomain  to  the  other,  as  dictated  by  the 
smoothness  of  the  solution  and  determined  by  the  Multi-Resolution  Analysis 
of  Section  4. 

Each  WENO  subdomain  must  have  a  ’’Buffer  Area”  in  order  to  detect  outgo¬ 
ing  shocks  or  high  gradients  and  ’’warn”  the  immediate  neighboring  spectral 
subdomains  to  switch  to  WENO  type.  The  ’’Buffer  Area”  for  a  given  WENO 
subdomain  is  dehned  as  two  sets  of  grid  points  at  either  ends  of  the  subdo¬ 
main,  namely,  Bq  =  {xq,  . . . ,  and  . . . ,  As  greater 

as  is  the  value  of  the  number  of  ghost  points,  Ng,  earlier  is  the  detection  of 
shocks  and  gradients,  however,  at  the  cost  of  earlier  switching  of  the  adjacent 
spectral  subdomains  and  greater  chance  of  unnecessary  costly  computations. 
Satisfactory  results  have  been  obtained  with  the  default  value  Ng  =  r,  the 
number  of  ghost  cells  used  in  the  given  order  of  the  WENO  scheme  in  this 
study.  It  should  be  noted  that  the  grid  points  in  the  ’’Buffer  Area”  are  part 
of  the  interior  grid  points  and  should  not  be  confused  with  the  WENO  ghost 
points. 

For  the  switching  algorithm,  the  three  conditions  below  are  the  main  rules  to 
be  followed: 

(1)  If  a  subdomain  contains  high  gradients,  then  switch  its  spatial  discretiza¬ 
tion  to  (or  keep  it  as)  WENO; 

(2)  If  high  gradients  are  present  in  the  ’’Buffer  Areas”  of  neighboring  sub- 


21 


domains,  then  switch  the  current  subdomain  to  (or  keep  it  as)  a  WENO 
subdomain; 

(3)  In  any  other  case,  switch  the  subdomain  to  (or  keep  it  as)  a  spectral 
subdomain; 

The  hrst  condition  above  avoids  the  Gibbs  phenomenon,  keeping  the  disconti¬ 
nuities  inside  WENO  subdomains.  The  second  condition  ensures  the  switching 
to  WENO  subdomain  in  order  to  allow  only  WENO-to-WENO  transmission  of 
high  gradients.  The  third  condition  improves  the  numerical  efficiency,  since  it 
ensures  that  smooth  parts  of  the  solution  will  always  be  contained  in  spectral 
subdomains. 


Multi- Resolution  analysis  of  the  solution  is  performed  at  the  beginning  of  every 
step  of  the  associated  temporal  scheme,  in  our  case,  a  third  order  Runge-Kutta 
scheme,  as  described  in  Section  6  below.  At  each  subdomain  k,  we  dehne  the 
smoothness  flag  variable,  Flag^,  at  each  grid  point  Xi  excluding  the  ghost 
points,  as. 


Flag" 


1.  Mi  I  ^  ^-MR 

0,  otherwise 


f  —  0, . . . ,  , 


(48) 


where  df  are  the  MR  coefficients.  Since  the  Multi- Resolution  Analysis  requires 
uniformly  spaced  grids,  the  spectral  grids  are  first  interpolated  to  uniformly 
spaced  grids  before  obtaining  the  MR  coefficients  df.  The  necessary  ghost 
points  are  acquired  from  the  neighboring  subdomains.  At  the  boundary  sub- 
domains,  the  values  of  the  boundary  ghost  points  are  extrapolated  linearly 
from  the  interior  data. 


The  algorithm  proceeds  by  checking  for  each  spectral  subdomain  k  and  at  the 
Buffer  Areas  of  the  neighboring  subdomains  k— 1  and  k-|-l,  if  any  of  {Flagf ,  i  = 
0, . . . ,  A^^},  {Flag^+\  f  =  0, . . . ,  Ng}  or  {Flagf"\  i  =  -  Ng  +  1, . . . ,  N^} 

is  equal  to  one.  If  so,  it  switches  subdomain  k  to  a  WENO  discretization. 
Otherwise,  a  spectral  discretization  is  implemented,  or  kept.  These  switches 
require  the  use  of  interpolation  from  a  Chebyshev  grid  to  a  uniformly  spaced 
one  and  vice-versa: 

•  To  switch  from  the  spectral  subdomain  to  the  WENO  subdomain,  the  data 
are  interpolated  onto  the  uniformly  spaced  grid  via  the  spectral  interpola¬ 
tion  formula. 

•  To  switch  from  the  WENO  subdomain  to  the  spectral  subdomain,  the 
data  are  interpolated  onto  the  Chebyshev  Gauss-Lobatto  points  via  the  La- 
grangian  interpolation  polynomial  of  the  same  order  as  the  WENO  method. 

Remark  3  Back  and  forth  switching  between  WENO  and  spectral  discretiza¬ 
tions  may  occur  too  frequently  for  the  same  domain  when  the  Cmr  is  marginally 
set.  The  d\  coefficients  might  oscillate  around  the  parameter  Cmr  in  time  due 
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to  some  numerical  factors  such  as  dissipation,  dispersion  and  nonlinear  ef¬ 
fects,  or  any  combination  of  such.  This  pattern  of  switching  can  repeat  itself 
for  a  while  until  the  solution  settles  down  with  a  clear  definition  of  the  df, 
which  is  either  greater  than  or  smaller  than  the  MR  tolerance  Cmr-  In  order  to 
alleviate  such  occurrences,  one  must  devise  a  procedure  preventing  the  switch 
from  WENO  to  spectral  if  it  had  already  occurred  recently.  However,  such  pro¬ 
cedure  must  never  prevent  a  spectral  to  WENO  switch,  for  oscillations  and 
instability  might  occur. 

We  end  this  section  proceeding  fnrther  with  the  temporal  integration  of  the 
inviscid  Bnrgers  Eqnation  (47)  to  a  final  time  t  =  5  where  the  shock  has  al¬ 
ready  developed  and  performed  an  entire  revolntion  at  the  periodical  domain. 
The  numerical  results  are  shown  at  figure  6  where  we  now  have  partitioned 
the  interval  with  10  subdomains.  Each  spectral  subdomain  uses  17  points, 
while  the  WENO  ones  use  40  points.  Note  that  the  Hybrid  Method  was  able 
to  compute  the  exact  location  of  the  shock,  demonstrating  its  conservative 
property. 


Fig.  6.  Inviscid  Burgers  Equation  solution;  The  solution  computed  by  the  Hybrid 
scheme  is  plotted  with  symbols  against  the  exact  solution  in  solid  line.  The  er¬ 
ror  is  plotted  in  dashed  lines.  Vertical  dashed  lines  show  the  subdomains  division. 
Subdomains  [0.2, 0.4],  [0.4, 0.6]  and  [0.6, 0.8]  are  WENO,  all  others  are  Spectral. 


6  Numerical  Experiments 


In  this  section  we  present  numerical  experiments  with  the  system  of  Euler 
Equations  for  gas  dynamics  in  strong  conservation  form: 

Qt  +  Fa,  =  0,  (49) 


where 
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Q  =  (p,  im,  Ef,  F  =  (p«,  pu^  +  P,  (F  +  P)uY , 


(50) 


and  the  equation  of  state 

P  =  -^){e +  \pu^)  ,  7  =  1.4.  (51) 


In  all  examples  below,  the  physical  domain  is  partitioned  into  a  fixed  number 
of  same  size  subdomains  and  the  initial  conhguration  of  these  subdomains 
depends  on  the  initial  condition  under  consideration.  Spectral  subdomains 
are  discretized  with  Chebyshev-Gauss-Lobatto  collocation  points  and  WENO 
subdomains  use  an  uniform  grid,  where  the  classical  fifth  order  characteristic- 
wise  WENO  finite  difference  is  applied.  We  shall  use  the  same  number  of 
uniformly  spaced  grid  points  for  all  WENO  subdomains,  as  well  as  the 
same  number  of  Chebyshev  collocation  points  Ng  at  all  spectral  subdomains, 
unless  noted  otherwise.  The  detection  of  discontinuities  and  high  gradients 
is  performed  through  a  hfth  order  Multi-Resolution  Analysis  applied  to  the 
density  function  p.  A  16_th  order  Exponential  hlter  is  employed  in  all  spectral 
subdomains. 

To  evolve  in  time  the  ODEs  resulting  from  the  semi-discretized  PDEs,  the 
third  order  Total  Variation  Diminishing  Runge-Kutta  scheme  (RK-TVD)  will 
be  used  [24]: 

U^  =  U'^  +  AtL{U'^) 

U^  =  ^(3U^ +  U^  +  AtL{U^))  ,  (52) 

f/^+i  =  +  2U^  +  2AtL{U^)) 

where  L  is  the  spatial  operator.  CEL  numbers  for  the  spectral  and  WENO 
subdomains  are  set  to  be  1  and  0.4,  respectively.  All  results  of  the  hybrid 
method  are  compared  with  an  exact  solution  or  with  a  numerical  solution 
obtained  with  the  classical  WENO  scheme  with  the  same  spatial  resolution  of 
the  hybrid  scheme.  This  is  achieved  by  using  the  same  number  of  points  for 
the  classical  scheme  as  if  all  subdomains  in  the  hybrid  were  of  the  WENO  type 
(see  Section  5).  The  main  goal  of  the  numerical  experiments  below  is  to  show 
that  the  Hybrid  method  obtains  equivalent  solutions  as  the  classical  WENO 
scheme,  however  at  a  lower  computational  cost.  This  is  easily  concluded  from 
the  facts  that  spectral  discretization  is  much  more  efficient  than  the  finite 
differences  at  smooth  parts  of  the  solution  and  that  spectral  subdomains  avoid 
the  expensive  characteristic  decomposition  of  the  WENO  scheme. 

In  the  figures  shown  for  the  one  dimensional  test  cases,  the  spectral  and 
WENO  subdomains  are  those  denoted  with  red  squares  (□)  and  black  trian¬ 
gles  (A)  symbols  respectively.  The  subdomains  are  also  indicated  by  vertical 
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dashed  lines  and  major  ticks  on  the  x  axis.  Under  each  hgnre,  the  snbdo- 
mains  conhgnration  is  indicated  as  a  seqnence  of  powers  of  S  and  W,  where 
the  exponent  means  the  nnmber  of  consecntive  snbdomains  of  the  same  type. 
This  notation  will  make  easier  the  connting  of  snbdomain  types  when  a  large 
nnmber  of  snbdomains  is  used. 


6.1  SOD  Tube  Problem 


We  hrst  consider  the  Sod  problem  of  the  Compressible  Euler  equations  with 
initial  Riemann  data: 


(P,  U,  P) 


(  0.125,  0,  0.1  )  -5  <  X  <  0 
(  1,  0,  1  )  0  <  X  <  5 


(53) 


The  physical  domain  [—5,  5]  is  partitioned  into  5  equal  subdomains.  Since  the 
initial  condition  is  discontinuous  at  the  center  of  the  physical  domain,  the 
middle  subdomain,  [—1, 1],  is  set  as  WENO  (Figure  7(a)).  Here,  Ng  =  16  and 
=  110  and  the  MR  tolerance  was  taken  as  =  5  x  10“^. 
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(a)  t  =  0  , 


(b)  t  =  0.39 


(d)  t  =  1.46 
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(c)  t  =  0.54 


(e)  t  =  2 


Fig.  7.  The  density  profile  of  the  Sod  shock  tube  problem  at  times  (a)  t  =  0,  (b) 
t  =  0.39,  (c)  t  =  0.54,  (d)  t  =  1.46  and  (e)  t  =  2  with  5  subdomains.  Ns  =  16, 
Nyy  =  110  and  Emr  =  5  x  10“^.  Spectral  (□);  WENO  (A);  The  solid  line  is  the 
exact  solution. 

As  time  evolves,  the  initial  density  discontinuity  develops  into  two  jump  dis¬ 
continuities,  the  leftward  moving  shock  front  and  contact  discontinuity  and  a 
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rightward  moving  rarefaction  wave,  which  is  discontinuous  in  the  hrst  deriva¬ 
tive  (Figures  7(b)  and  (c)).  Notice  that  when  the  discontinuities  move  on  closer 
to  the  boundary  of  the  WENO  subdomain  [—1,1],  the  neighboring  spectral 
subdomains  are  switched  to  WENO  ones.  At  a  later  time,  the  discontinuities 
have  moved  further  to  the  left  side  of  the  physical  domain,  however  the  hrst 
and  last  subdomains  were  kept  as  spectral  ones,  since  none  of  them  has  been 
’’threatened”  by  the  discontinuities. 


6.2  123  Tube  Problem 


We  now  consider  the  123  problem  with  initial  Riemann  data. 


(  1,  -2,  0.4  )  -5  <  a:  <  0 
ip,U,P)  =  {  ^  ^  - 

(1,  2,0.4)  0<a:<5 


(54) 


The  solution  consists  of  two  rarefaction  waves  moving  in  opposite  directions 
which  are  generated  at  the  center  of  the  physical  domain  by  a  discontinuity 
in  the  velocity. 


(a)  f  =  0  s^w^si 


(b)  t  =  0.25 


(c)  t  =  0.29 
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(d)  t  =  0.62  s^wisiwis^  (e)  t  =  1  siwis^wisi 


Fig.  8.  The  density  profile  of  the  123  problem  at  times  (a)  t  =  0,  (b)  t  =  0.25,  (c) 
t  =  0.29,  (d)  t  =  0.62  and  (e)  t  =  1  with  5  subdomains.  Ng  =  16,  =  110  and 

Cma  =  5  X  10“^.  Spectral  (□);  WENO  (A);  The  solid  line  is  the  exact  solution. 


Using  the  same  setting  as  the  Sod  problem  discussed  in  the  previous  example, 
we  compute  the  numerical  solution  up  to  t  =  1.  Figure  8(d)  shows  that  the 
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middle  WENO  subdomain  [-1,1]  is  switched  to  spectral  as  soon  as  the  rar¬ 
efaction  waves  move  away  from  the  buffer  zones  of  the  neighboring  domains. 
Note  also  that  the  hrst  and  last  spectral  subdomains  have  a  smaller  number 
of  discretization  points  than  the  newly  created  spectral  subdomain,  showing 
that  the  Hybrid  scheme  might  also  provide  quantitative  flexibility  at  the  lo¬ 
cal  discretizations.  Non-consecutive  WENO  subdomains  can  also  have  distinct 
number  of  grid  points. 


6.3  Blastwaves  Simulation 


The  one  dimensional  Blast  waves  interaction  problem  by  Woodward  and  Col- 
lela  [17]  has  the  following  initial  prohle 


(1,0, 

1000 ; 

)  0  <  x  <  0.1 

{p,u,p)  =  < 

(1,0, 

0.01 ; 

)  0.1  <  a;  <  0.9  • 

(55) 

(1,0, 

100 ; 

)  0.9  <  a;  <  1.0 

The  initial  pressure  gradients  generate  two  density  shock  waves  that  collide 
and  interact  later  in  time.  In  this  experiment,  the  physical  domain  is  subdi¬ 
vided  into  10  subdomains  and  the  number  of  Chebyshev  collocation  points 
and  WENO  uniform  grid  are  Ng  =  16  and  =  50,  respectively.  Reflective 
boundary  conditions  are  applied  at  both  ends  of  the  physical  domain.  Figure  9 
shows  that  at  the  initial  times,  the  expensive  WENO  discretizations  are  local¬ 
ized  around  the  density  peaks,  while  the  remaining  subdomains  are  spectral. 
At  intermediate  times,  the  left  rarefaction  wave  spreads  around,  requiring  the 
use  of  more  WENO  subdomains,  however,  by  the  end  of  the  simulation,  we 
have  more  spectral  subdomains  than  WENO  ones.  It  is  clear  from  hgure  9(c) 
that  one  could  decrease  the  number  of  discretization  points  at  the  hrst  hve 
spectral  subdomains,  or  take  a  step  further  and  merge  all  the  subdomains 
at  a  single  spectral  one,  increasing  the  numerical  efficiency  of  the  algorithm. 
These  improvements  will  be  considered  in  future  implementations  of  the  Hy¬ 
brid  method. 


27 


(a)  t  =  0.014  (b)  t  =  0.16  (c)  t  =  0.38 


Fig.  9.  The  density  profile  of  the  Blastwave  problem  at  times  t  =  0.014,  t  =  0.16 
and  t  =  0.38  with  10  subdomains  as  indicated  by  the  vertical  dash  lines,  A^s  =  16, 
Nw  =  50  and  Emr  =  5  x  10“^.  Spectral  (□);  WENO  (A);  The  solid  line  is  the 
solution  computed  using  the  classical  WENO  scheme  with  500  points. 

6-4  Shock- Entropy  Wave  Interaction 


Consider  the  one  dimensional  Mach  3  shock-entropy  wave  interaction,  specihed 
by  the  following  initial  conditions: 

f  (  3.857143,  2.629369,  10.33333  )  -5  <  x  <  -4 

(p,n,P)  =  r  '  -  ,  (56) 

[^  (  1  +  esin(A:a;),  0,  1  )  — 4  <  a;  <  15 

where  x  G  [—5, 15] ,  e  =  0.2  and  k  =  5.  The  solution  of  this  problem  consists 
of  shocklets  and  hne  scales  structures  which  are  located  behind  a  right-going 
main  shock.  Figure  10  shows  that  the  hybrid  method  is  able  to  capture  all 
the  smooth  high  frequency  waves  behind  the  main  shock  with  spectral  dis¬ 
cretizations.  Note  that  the  WENO  subdomains  are  located  only  at  the  main 
shock  and  at  the  steep  gradients  of  the  N-waves.  The  hybrid  method  uses  40 
subdomains,  with  =  16  and  Ny^r  =  50.  The  solid  black  hne  is  the  solution 
computed  with  the  hfth  order  WENO  scheme  with  2000  grid  points.  Figures 
10  (d)  also  shows  that  even  with  the  great  complexity  of  the  solution,  less 
than  30%  of  the  total  number  of  subdomains  are  of  the  WENO  type. 


6. 5  Shock-  Vortex  Interaction 


We  now  show  some  previous  results  of  the  Hybrid  method  in  more  spatial 
dimensions.  A  more  detailed  discussion  will  be  reported  at  an  upcoming  arti¬ 
cle.  In  this  experiment,  the  multidomain  hybrid  conhguration  is  applied  to  a 
2D  Shock- Vortex  interaction.  We  consider  a  counter-clockwise  rotating  vortex 
centered  at  (xc,  Vc),  and  strength  F,  with  a  tangential  velocity  prohle  [8]  given 
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Fig.  10.  The  density  profile  of  the  Mach  3  Shock-Entropy  wave  interaction  at  times 
(a)  t  =  0,  (b)  t  =  1.5,  (c)  t  =  3.2  and  (d)  t  =  4.5  with  40  subdomains,  as  indicated 
by  the  vertical  lines.  Ng  =  16,  =  50  and  Cmh  =  5  x  10“^.  Spectral  (□);  WENO 

(A);  The  solid  line  is  the  solution  computed  with  the  hfth  order  WENO  scheme 
with  2000  grid  points. 

in  polar  coordinates  by: 


rr(ro  ^  0  <  r  <  ro  <  ri 


f/(r)  =  < 


rr(r  ^ 


Ti  ro  <  T  <  ri 


0 


r  >  n 


(57) 


where  tq  =  0.2  and  ri  =  1.0.  The  physical  domain  (0  <  a;  <  3.9,  —2  <  y  <  2)  is 
partitioned  into  a  13  x  10  grid  of  subdomains.  Spectral  and  WENO  subdomains 
use  a  32  X  32  grid  of  Chebyshev  collocation  points  and  a  50  x  50  uniform 
grid,  respectively.  The  Multi-Resolution  analysis  is  performed  with  tolerance 


We  simulate  a  Mach  3  shock  interacting  with  the  vortex  of  strength  T  =  0.25 
and  compare  the  solution  obtained  by  the  hybrid  method  with  the  above  setup 
with  a  highly  resolved  solution  using  the  classical  fifth  order  characteristic- wise 
WENO  hnite  difference  scheme  with  a  1200  x  1200  grid  at  hnal  time  t  =  0.6 
(Figure  11).  The  shock  and  the  high  gradient  regions  immediately  behind 
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the  shock  are  well  captured  by  the  MR  analysis,  as  indicated  by  the  WENO 
subdomains  enclosed  with  black  bounding  boxes.  The  remaining  subdomains 
are  accurately  dealt  with  by  the  spectral  method  since  all  the  essential  small 
and  large  scale  structures  of  the  flow  held  are  correctly  represented  in  the 
hybrid  solution.  Once  again,  the  number  of  WENO  subdomains  is  far  fewer 
than  the  spectral  ones,  resulting  in  a  more  efficient  algorithm  than  the  classical 
WENO  scheme. 


Fig.  11.  Density  contour  of  the  Shock- Vortex  interaction  with  Mach  number  Mg  =  3 
and  vortex  strength  F  =  0.25  at  t  =  0.6.  Hybrid  (Left);  Classical  WENO  (Right). 
WENO  subdomains  are  denoted  by  bounding  boxes. 


7  Conclusions 


We  have  presented  the  one  dimensional  version  of  the  multi-domain  spatial 
and  temporal  adaptive  Hybrid  Spectral- WENO  method  for  nonlinear  systems 
of  hyperbolic  conservation  laws.  The  main  idea  is  to  partition  the  physical 
domain  into  a  grid  of  subdomains  which  are  either  a  Chebyshev  collocation 
grid  for  the  spectral  method  or  an  uniformly  spaced  grid  for  the  high  order 
WENO  method.  High  order  Multi- Resolution  Analysis  is  used  to  measure 
the  smoothness  of  the  solution  in  a  given  subdomain  in  order  to  switch  a 
subdomain  from  spectral  to  WENO  when  the  solution  becomes  non-smooth, 
and  vice  versa.  In  this  work,  the  Hybrid  method  was  tested  with  the  standard 
shock-tube  test  problems,  Shock-Entropy  wave  interaction  problem  and  the 
Blastwave  problem.  The  results  matched  the  ones  computed  by  the  classical 
hfth  order  WENO  hnite  difference  method.  Timing  results  will  be  given  in 
the  upcoming  paper  in  the  two  dimensional  extension  of  the  Hybrid  method, 
which  is  a  much  more  meaningful  measure  than  the  one  dimensional  one. 

Currently,  we  are  investigating  the  role  of  the  Multi-Resolution  tolerance 
parametereMfl  and  its  effects  on  the  Hybrid  solution.  Furthermore,  the  imple¬ 
mentation  of  the  smoothness  measurements  in  the  spectral  subdomains  can 
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be  improved  and  is  under  investigation.  We  plan  to  extend  the  Hybrid  method 
to  nonlinear  hyperbolic  conservation  laws  system  in  higher  dimensions  with 
higher  orders  WENO  schemes. 
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